Excited stationary states of trapped Bose-Einstein condensates 
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We investigate the excited stationary states of Bose-Einstein condensates trapped in harmonic 
potentials. We derive simple analytical approximations of the first few eigenstates of the associated 
time-independent one-dimensional Gross-Pitaevskii equation and their energies. Our results are 
excited state generalizations of the Thomas-Fermi approximation of the ground state. 
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As opposed to traditional studies of liquid helium, the 
recent experimental work on Bose-Einstein condensation 
centers on inhomogeneous mesoscopic systems. How- 
ever, such systems are highly nonlinear and no exact an- 
alytical expressions for the wavefunctions and associate 
energies are yet known. In order to understand the be- 
havior of Bose-Einstein condensates (BEC) it would be 
beneficial to have at least simple approximate expressions 
for the condensate wave functions. These could be used 
to study the qualitative behavior of the condensate and, 
given the approximate expressions are precise enough, 
might even be useful for the determination of overlap fac- 
tors between different wave functions, spectra and other 
quantities of interest. Yukalov et al. ^ have investigated 
this problem using a perturbation method. Also, Kivshar 
et al. ^ recently highlighted the connection between this 
problem and that of dark solitons and gave numerical 
solutions to the first few excited stationary states. In 
this communication we derive simple analytical approx- 
imations to the stationary states $„ (n = 0, 1,2, . . .) of 
the one-dimensional Gross-Pitaevskii equation for a BEC 
confined in a harmonic trap |4| : 
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Here z is the longitudinal trap coordinate, the aspect 
ratio, A = uJz/^ ^ 1, is given by the ratio of the lon- 
gitudinal trap frequency ojz to the transversal trap fre- 
quency uj, A formally is the effective spring constant of 
the 1-dimensional problem. A is assumed to be much 
smaller than unity in order to yield a cigar-shaped con- 
densate which can be modeled by a 1-dimensional prob- 
lem. Q/2 — AnaN/ao is the effective non-linearity for 
the interaction of the N trapped atoms with mass m and 
scattering length a in the trap with an effective diameter 

The states $„ are normalized to unity, but, owing to 
the non-linearity of the Gross-Pitaevskii-equation, they 
do not form an orthogonal set, generally ($„|<i>m) ^ 0. 

It is convenient to rescale the Gross-Pitaevskii- 
equation to make the coefficients of the kinetic and the 
non-linear term equal to unity. Using the rescaling 
Yn = <&„ i -^/Q/eo and x = Zs/Ieq Eq. (|l|) becomes 



The advantage of this scaling is that the solutions are 
parameterized solely by A. The energy eigenvalue e„ — 
Sn/so is now scaled in units of the Thomas- Fermi en- 
ergy £o = (3^/^/8) • [QXy/^ of the lowest stationary 
state and the spring constant is expressed in terms of 
A^ = (wz/w)^ = (512/9) • Eq/Q^. Our choice of rescaling 



results in the normalization ($n|<&„) — 1 
which can be written as 
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The Thomas-Fermi approximation of Yq entirely ne- 
glects the contribution from the kinetic energy term and 
yields the standard approximate solution M of (0): 



Y^{x) = Vr^A2lE2^ for \x\ <Xo = 1/A (4) 
and Yo{x) — otherwise, as illustrated in Fig. |l]. 




FIG. 1. Plot of the parabolic trap potential, the standard 
Thomas-Fermi solution Q for the condensate state Yo (dashed 
line) and the corresponding exact numerical solution (solid 
line). The value of the nonlinearity is A = 1/30. Both solu- 
tions have been displaced vertically by their respective energy 
values and share the same arbitrary vertical scaling. 
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To generalize the Thomas-Fermi approach to excited 
stationary states, F„ for n > 0, we need to include the 
kinetic energy term. Our approach to deriving a solution 
is based on a piecewise ansatz which is motivated as fol- 
lows. Over small regions about some fixed point x — x' 
the potential energy term X^x^ can be considered slowly 
varying compared to the other terms in Eq. (^. The lo- 
cal solution about x' is then given approximately by the 
solution of the nonlinear Schrodinger equation 



ox^ 



Ynix) - 



(5) 



where = e„ — x'^ . This equation is well known 
in soliton theory; the solutions with finite amplitude are 
given by the Jacobi-sinc function 'sn': 



A sn 



(^^e'„{l^Ay2){x - a), 



A 



(6) 



where a is an arbitrary displacement along the x axis and 
< A < 1. The wavelength of the function in is given 
by 



A = 



K{ 



) 



(7) 



where iiT is a complete elliptic integral of the first kind 
[ pT| . Imagine joining the solutions for two different val- 
ues of x' smoothly at some intermediate point along the x 
axis. As an example, let this intermediate point be a zero 
of the two solutions. The slope of (^) at a zero is given by 
e^Aa/I — A'^/2 and so a smooth join requires the value of 
A to be larger for the solution with the larger value of Ix'l 
to compensate for the smaller value of = e„ — X^x'^. 
Repeating this smooth joining of solutions for successive 
pairs of x' values along the a: axis leads to a set of piece- 
wise solutions for which the value of A increases with 
\x'\. At the critical value A = 1 the wavelength becomes 
infinite and the solution of Eq. (^) is given by the corre- 
sponding limit of m) namely 
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This piecewise analysis implies that the general finite- 
amplitude solution of Eq. (||) will be oscillatory near 
the origin and change to a tanh-like function some dis- 
tance away. Moreover, the tanh function in (H) is rel- 
atively flat for \x — a\ significantly different from zero, 
and so (p|) is approximated well by just the prefactor 
^/e^^ — v^n — A^x'^. We note that this prefactor is the 
Thomas-Fermi solution (^) in terms of the local variable 
x'. 

We can reduce the number of piecewise segments 
needed in our ansatz using energy considerations ^ . The 
potential energy and its variation are smallest in the cen- 
ter of the trap and to a good approximation we can ne- 
glect the trap potential in this region. This amounts to 



replacing with e„ and holding A constant in our os- 
cillatory solution (H) to give an oscillating function Yosc 
which has a fixed wavelength and amplitude. Further 
away from the center of the trap the potential energy 
increases at the expense of the kinetic energy to the ex- 
tent that the later becomes negligible. The wave function 
is then well described by the conventional Thomas-Fermi 
solution \/en — X^x!^ . Thus, replacing x' with x in the 
prefactor Ve„ — X^x''^ of (||) yields a function Ytf which 
provides a transition between Yosc and the Thomas-Fermi 
solution. We will call Ytf the "tail function" as it gives 
the shape of the condensate in its outer regions. 

Our ansatz is to match these two segments l^sc and 
Itf smoothly at points x — ±T„ which lie symmetri- 
cally about the origin. In analogy to Eq. we also 
set Yn[x) — for \x\ > X = ^fl^jX where X is the 
half-width of the condensate in the Thomas-Fermi ap- 
proximation. Thus our ansatz is given by 
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where 
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The points x — ±T„, where the solution switches be- 
tween the oscillatory part and the tail function is given 
by r„ = (n - 1) A/4. The forms of Tn and Y^scix) en- 
sure that Yn{x) is either symmetric or antisymmetric for 
n even or odd. 

The parameters A, T„ and e„ are determined by the 
requirements that the two segments match smoothly and 
that the solution is normalized according to Eq. 
The first requirement is satisfied when the gradients, 

T„ are equal, i.e. when 



i^Yoso\x=T„ and ^Ftf 
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For notational simplicity we concentrate on the second 
excited state I2, i-e. n — 2. The expression for T2 can be 
approximated by T2 « ■y/2/e2 arctanh(-\/ (1 + A^)/2) « 
\/2/e2-ln(2/Vl - A). This can be used to solve Eq. ^ 
for A « 1 - A £2 VF(8e2/A)/2 where W is the Lambert 
function defined by W^exp(W^) = x. This, in turn, can 
be approximated to give A k, \ ~ X £2 ln(8e2/A)/2 and 
therefore T2 « - ln(A £2 ln(8e2/A)/8)/V5^- Thus we 
now have a simple expression for T2 in terms of £2 . 



2 



The second requirement, i.e. that Y2{x) be normalized 
according to Eq. (||), determines the value of 62- The con- 
tribution from the oscillatory part of the wave function 



(^oscli^osc) ^2 V2 arctanh(v/(l+A2)/2) 



£2 

and the contribution from the tails is approximately 
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+ (e„ - A2T2)tanh2(V(e„ - X^x^)/2)\ . (11) 

We use tanh(-^0£) w 1 to find a simplified expression 
for (p]). Inserting the above expression for A and T 
into (1^21^2) = 4/(3A), treating A and 62 ~ £2 ~ eo as 
small expansion parameters and performing a somewhat 
tedious calculation using an expansion in 62 followed by 
an expansion in A finally yields 62 ~ 1 + 2^2 A. 

The derivation for the general case n — 1,2,3,... is 
very similar. We find that 



and 



e„ « 1 + nV2 A , 
^« 1- Ae„ln(8e„/A)/2 , 
T„ « (1 - n) ln(A e„ ln(8e„/A)/8)/> 



(12) 
(13) 
(14) 



Eq. (|) together with Eqs. are the main results of 

this communication and give our analytical approxima- 
tions of the stationary states of a harmonically trapped 
condensate. We note, however, that the greater the value 
of n, the more stringent the requirements on the aspect 
ratio A to be small. 

We should expect, and do obtain, best results for low 
order n because in this case the above approximations 
are least critical. In a more detailed paper [|| we will 
present a self-consistency argument to justify our ap- 
proach. Figs. ||-^ compare our analytical and the exact 
(numerically determined) wave functions for low-order 
excited states and illustrate the increasing deviation from 
the fixed wavelength approximation for more highly ex- 
cited states. As in the Thomas-Fermi solution in Fig. ||, 
the omission of the kinetic energy for the outer regions 
of the condensate leads to the artificial truncation at the 
outer edge x — X of the wave functions rather than a 
Gaussian tail. Also Fig. ^ shows that the fixed wave- 
length approximation of Eq. (^) for the oscillatory seg- 
ment of the solution breaks down for high excitation 
numbers. Interestingly we have found that the overlap 
between the numerical and the analytical solutions to I4 
is higher for larger values of A (e.g. A = 1/20). We 
believe this is a consequence of the fixed wavelength ap- 
proximation, we will explore this point elsewhere. 




FIG. 2. Plot of the parabolic trap potential, our ap- 
proximate analytical solution for the condensate state Yi{x) 
(dashed line) and the corresponding numerical solution (solid 
line) for a nonlinearity parameter of A = 1/30. As in Fig. ^ 
both solutions have been displaced vertically by their respec- 
tive energy values and they share the same arbitrary vertical 
scaling. 




FIG. 3. Analogous to Fig. ^ for the second excited station- 
ary state Y2{x). 
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FIG. 4. Analogous to Fig. ^ for the fourth excited station- 
ary state Y4,{x). 

These excited stationary states may soon be realized 
experimentally. Various possibilities to generate excited 
states have been discussed in the literature ■ More- 

over, Burger et al. ^ recently produced a slowly-moving 
dark soliton using the 'phase imprinting method' [^,^. 
It should be possible with a slight modification of their 
experiment to generate the first excited state Yi{x) as a 
stationary dark soliton. The generation of higher excited 
states should also be possible in a similar manner. 

In conclusion, wc have derived simple, approximate 
analytic expressions of the excited stationary solutions 
of the 1-dimensional Gross-Pitaevskii equation using an 
ansatz based on piecewise solutions. Our analytic solu- 
tions agree very well with numerical solutions for a suf- 
ficiently small aspect ratio. Generalizations of the ap- 
proach presented here should be rather straightforward, 
for example, one can introduce a coordinate dependent 
wavelength of the oscillatory part of the solution rather 



than using the fixed wavelength approximation. Further 
details will be explored elsewhere 
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